00001 /** Functions for SCS inverse and division 00002 00003 @file division_scs.c 00004 00005 @author Defour David David.Defour@ens-lyon.fr 00006 @author Florent de Dinechin Florent.de.Dinechin@ens-lyon.fr 00007 00008 This file is part of the SCS library. 00009 00010 00011 */ 00012 00013 /* 00014 Copyright (C) 2002 David Defour and Florent de Dinechin 00015 00016 This library is free software; you can redistribute it and/or 00017 modify it under the terms of the GNU Lesser General Public 00018 License as published by the Free Software Foundation; either 00019 version 2.1 of the License, or (at your option) any later version. 00020 00021 This library is distributed in the hope that it will be useful, 00022 but WITHOUT ANY WARRANTY; without even the implied warranty of 00023 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00024 Lesser General Public License for more details. 00025 00026 You should have received a copy of the GNU Lesser General Public 00027 License along with this library; if not, write to the Free Software 00028 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00029 00030 */ 00031 #include "scs.h" 00032 #include "scs_private.h" 00033 00034 00035 /* 00036 * Compute 1/x with a Newton scheme 00037 */ 00038 void scs_inv(scs_ptr result, scs_ptr x){ 00039 scs_t tmp, res, res1, scstwo; 00040 double app_x, inv; 00041 00042 scs_set(tmp, x); tmp->index = 0; scs_get_d(&app_x, tmp); 00043 00044 scs_set_si(scstwo, 2); 00045 /* inv is a 53-bit approximation of 1/x */ 00046 inv = 1/app_x; 00047 00048 scs_set_d(res, inv); 00049 res->index -= x->index; 00050 00051 /* First Newton Iteration */ 00052 scs_mul(res1, x, res); 00053 scs_sub(res1, scstwo, res1); 00054 scs_mul(res, res, res1); 00055 00056 /* Second Newton Iteration */ 00057 scs_mul(res1, x, res); 00058 scs_sub(res1, scstwo, res1); 00059 scs_mul(result, res, res1); 00060 00061 return; 00062 } 00063 00064 /* 00065 * Compute result = x/y; 00066 */ 00067 void scs_div(scs_ptr result, scs_ptr x, scs_ptr y){ 00068 scs_t res; 00069 00070 if (X_EXP != 1){ 00071 R_EXP = X_EXP / Y_EXP; 00072 return; 00073 } 00074 00075 scs_inv(res, y); 00076 scs_mul(result, res, x); 00077 return; 00078 }